﻿// SyntheticDataGenerator.Generators.PoissonGenerator
//
// (c) 2011 Arthur Pitman
//
// Part of Market-Basket Synthetic Data Generator
// A C# adaptation of the IBM Quest Market-Basket Synthetic Data Generator
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// (Version 2) as published by the Free Software Foundation.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public
// License (Version 2) along with this program; if not, write to 
// the Free Software Foundation, Inc., 
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
// or see <http://www.gnu.org/licenses/>.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SyntheticDataGenerator.Generators
{
	/// <summary>
	/// Generates random numbers with a poisson distribution
	/// 
	/// See "Numerical Recipes in C: The Art of Scientific Computing" by Press et al.
	/// </summary>
	public class PoissonGenerator : IGenerator
	{
		/// <summary>
		/// Mean and variance of the distribution
		/// </summary>
		protected double lambda;

		/// <summary>
		/// Algorithm parameters
		/// </summary>
		protected double sq, alxm, g;

		/// <summary>
		/// The <see cref="UniformGenerator"/> used for randomness
		/// </summary>
		protected UniformGenerator uniformGenerator = new UniformGenerator();

		/// <summary>
		/// Algorithm coefficients
		/// </summary>
		protected static double[] coefficients= new double[]
		{
			76.18009173,
			-86.50532033,
			24.01409822,
			-1.231739516,
			0.120858003e-2,
			-0.536382e-5
		};

		/// <summary>
		/// Initializes a new <see cref="PoissonGenerator"/> instance with the specified mean and variance
		/// </summary>
		/// <param name="lambda">The mean and variance of the poisson distribution</param>
		public PoissonGenerator(double lambda)
		{
			this.lambda = lambda;        

			// initialize algorithm parameters
			if (lambda < 12.0)
			{
				g = Math.Exp(-lambda);
			}
			else
			{
				sq = Math.Sqrt(2.0 * lambda);
				alxm = Math.Log(lambda);
				g = lambda * alxm - GammaLn(lambda + 1.0);
			}
		}

		/// <summary>
		/// Returns the next poisson distributed random value as an int
		/// </summary>
		/// <returns>random value</returns>
		public int NextInt()
		{
			return (int)Math.Round(Next());
		}

		/// <summary>
		/// Returns the next poisson distributed random value 
		/// </summary>
		/// <returns>random value</returns>
		public double Next()
		{
			double value;
			if (lambda < 12.0)
			{
				value = -1;
				double t = 1.0;

				do
				{
					value += 1.0;
					t *= uniformGenerator.Next();
				} 
				while (t > g);
			}
			else
			{
				double t, y;
				do
				{
					do
					{
						y = Math.Tan(Math.PI * uniformGenerator.Next());
						value = sq * y + lambda;
					} 
					while (value < 0.0);
					
					value = Math.Floor(value);
					t = 0.9 * (1.0 + y * y) * Math.Exp(value * alxm - GammaLn(value + 1.0) - g);
				} 
				while (uniformGenerator.Next() > t);
			}
			return value;
		}
	
		/// <summary>
		/// Calculates the logarithm of gamma function		
		/// See "Numerical Recipes in C" for more details of the algorithm
		/// </summary>
		/// <param name="x"></param>
		/// <returns></returns>
		protected double GammaLn(double x)
		{					
			double y = x - 1.0;			
			double z = y + 5.5;
			z -= (y + 0.5) * Math.Log(z);

			double ser = 1.0;
			for (int index = 0; index < 6; index++)
			{
				y += 1.0;
				ser += coefficients[index] / y;
			}

			return -z + Math.Log(2.50662827465 * ser);
		}
	}
}
